#define Sqrt sqrt
#define Power pow


/*
 *********************************************************
 ******************** Structures *************************
 *********************************************************
*/

struct parameters {
	double sdelta;
	double J;
	double t;
	double U;
	double beta;
	unsigned int L;
};

struct variables {
	double BDelta;
	double chi;
	double m;
	double d;
	double mu;
	double l_chi;
	double l_n;
	double l_m;
	double l_BDelta;
};

/*
 *********************************************************
 ***************** Small useful functions ****************
 *********************************************************
*/

void clear(parameters &p) {
	p.sdelta = p.J = p.t = p.J = p.beta = 0.;
	p.L = 0;
}

void clear(variables &v) {
	v.BDelta = v.chi = v.m = v.d = v.mu = 0;
	v.l_BDelta = v.l_chi = v.l_m = v.l_n = 0;
}

void addTo(variables &v, variables &q, double d=1.) {
	q.BDelta += v.BDelta*d;
	q.chi += v.chi*d;
	q.m += v.m*d;
	q.d += v.d*d;
	q.mu += v.mu*d;
	q.l_BDelta += v.l_BDelta*d;
	q.l_chi += v.l_chi*d;
	q.l_m += v.l_m*d;
	q.l_n += v.l_n*d;
}

void copyTo(variables &v, variables &q, double d=1.) {
	q.BDelta = v.BDelta*d;
	q.chi = v.chi*d;
	q.m = v.m*d;
	q.d = v.d*d;
	q.mu = v.mu*d;
	q.l_BDelta = v.l_BDelta*d;
	q.l_chi = v.l_chi*d;
	q.l_m = v.l_m*d;
	q.l_n = v.l_n*d;
}

inline double abs(double a) {
	if(a<0) a=-a;
	return a;
}

/*
 *********************************************************
 ***************** CLASS Functions ***********************
 *********************************************************
*/

class Func {
private:
	parameters p; // parameters
	variables v; // variables
	double n, r, w; // contractions

	// Abbreviations
	double BDelta;
	double chi;
	double d;
	double m;
	double mu;
	double lambdachi;
	double lambdan;
	double lambdam;
	double lambdaBDelta;
	double J;
	double t;
	double abgt;
	double abgs;

	void refresh() {
		n = 1.-p.sdelta;
		r = n/2. + v.m;
		w = n/2. - v.m;
		J = p.J;
		BDelta = v.BDelta;
		chi = v.chi;
		d = v.d;
		m = v.m;
		mu = v.mu;
		lambdachi = v.l_chi;
		lambdan = v.l_n;
		lambdam = v.l_m;
		lambdaBDelta = v.l_BDelta;
		J = p.J;
		t = p.t;
		abgt = gt();
		abgs = gs();
	}

public:
	Func(const parameters &p, const variables &v) : p(p), v(v) { refresh(); }

	void setPar(const parameters &p) { this->p = p; refresh(); }
	void setVar(const variables &v) { this->v = v; refresh(); }

	// Gutzwiller
	double gt();
	double gs();
	double dgtdm();
	double dgtdd();
	double dgtdn();
	double dgsdm();
	double dgsdd();
	double dgsdn();

	double gamma(double kx, double ky);
	double eta(double kx, double ky);
	double epsilon(double kx, double ky);
	double BDelta_d();
	double BDelta_af();

	// Energy
	double E1(double kx, double ky);
	double E3(double kx, double ky);

	double dE1dBDelta(double kx, double ky);
	double dE1dchi(double kx, double ky);
	double dE1dm(double kx, double ky);
	double dE1dgt(double kx, double ky);
	double dE1dgs(double kx, double ky);
	double dE1dmu(double kx, double ky);
	double dE1dl_chi(double kx, double ky);
	double dE1dl_n(double kx, double ky);
	double dE1dl_m(double kx, double ky);
	double dE1dl_BDelta(double kx, double ky);

	double dE3dBDelta(double kx, double ky);
	double dE3dchi(double kx, double ky);
	double dE3dm(double kx, double ky);
	double dE3dgt(double kx, double ky);
	double dE3dgs(double kx, double ky);
	double dE3dmu(double kx, double ky);
	double dE3dl_chi(double kx, double ky);
	double dE3dl_n(double kx, double ky);
	double dE3dl_m(double kx, double ky);
	double dE3dl_BDelta(double kx, double ky);
};

double Func::gt() {
	double part1 = (n-2.*v.d)/(n-2.*r*w);
	double part2 = sqrt((1.-w)*(1.-n+v.d)/(1.-r)) + sqrt(w*v.d/r);
	double part3 = sqrt((1.-r)*(1.-n+v.d)/(1.-w)) + sqrt(r*v.d/w);
	return part1*part2*part3;
}

double Func::gs() {
	return pow((n-2.*d)/(n-2.*w*r), 2.);
}

double Func::dgsdm() {
	return (-2*(-2*(-m + n/2.) + 2*(m + n/2.))*Power(-2*d + n,2))/Power(-2*(-m + n/2.)*(m + n/2.) + n,3);
}

double Func::dgsdd() {
	return (-4*(-2*d + n))/Power(-2*(-m + n/2.)*(m + n/2.) + n,2);
}

double Func::dgsdn() {
	return (-2*(1 - n)*Power(-2*d + n,2))/Power(-2*(-m + n/2.)*(m + n/2.) + n,3) +
			   (2*(-2*d + n))/Power(-2*(-m + n/2.)*(m + n/2.) + n,2);
}

double Func::dgtdm() {
	return -(((Sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.)) + Sqrt((d*(-m + n/2.))/(m + n/2.)))*
	        (-2*(-m + n/2.) + 2*(m + n/2.))*(Sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.)) +
	          Sqrt((d*(m + n/2.))/(-m + n/2.)))*(-2*d + n))/Power(-2*(-m + n/2.)*(m + n/2.) + n,2)) +
	   ((Sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.)) + Sqrt((d*(-m + n/2.))/(m + n/2.)))*
	      ((-(((1 + d - n)*(1 - m - n/2.))/Power(1 + m - n/2.,2)) - (1 + d - n)/(1 + m - n/2.))/
	         (2.*Sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.))) +
	        (d/(-m + n/2.) + (d*(m + n/2.))/Power(-m + n/2.,2))/(2.*Sqrt((d*(m + n/2.))/(-m + n/2.))))*(-2*d + n))/
	    (-2*(-m + n/2.)*(m + n/2.) + n) + ((((1 + d - n)/(1 - m - n/2.) +
	           ((1 + d - n)*(1 + m - n/2.))/Power(1 - m - n/2.,2))/
	         (2.*Sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.))) +
	        (-((d*(-m + n/2.))/Power(m + n/2.,2)) - d/(m + n/2.))/(2.*Sqrt((d*(-m + n/2.))/(m + n/2.))))*
	      (Sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.)) + Sqrt((d*(m + n/2.))/(-m + n/2.)))*(-2*d + n))/
	    (-2*(-m + n/2.)*(m + n/2.) + n);
}

double Func::dgtdd() {
	return (-2*(Sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.)) + Sqrt((d*(-m + n/2.))/(m + n/2.)))*
		      (Sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.)) + Sqrt((d*(m + n/2.))/(-m + n/2.))))/
		    (-2*(-m + n/2.)*(m + n/2.) + n) + ((Sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.)) +
		        Sqrt((d*(-m + n/2.))/(m + n/2.)))*((1 - m - n/2.)/
		         (2.*Sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.))*(1 + m - n/2.)) +
		        (m + n/2.)/(2.*(-m + n/2.)*Sqrt((d*(m + n/2.))/(-m + n/2.))))*(-2*d + n))/
		    (-2*(-m + n/2.)*(m + n/2.) + n) + (((1 + m - n/2.)/
		         (2.*(1 - m - n/2.)*Sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.))) +
		        (-m + n/2.)/(2.*Sqrt((d*(-m + n/2.))/(m + n/2.))*(m + n/2.)))*
		      (Sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.)) + Sqrt((d*(m + n/2.))/(-m + n/2.)))*(-2*d + n))/
		    (-2*(-m + n/2.)*(m + n/2.) + n);
}

double Func::dgtdn() {
	return -(((Sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.)) + Sqrt((d*(-m + n/2.))/(m + n/2.)))*
	        (Sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.)) + Sqrt((d*(m + n/2.))/(-m + n/2.)))*(1 - n)*
	        (-2*d + n))/Power(-2*(-m + n/2.)*(m + n/2.) + n,2)) +
	   ((Sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.)) + Sqrt((d*(-m + n/2.))/(m + n/2.)))*
	      (Sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.)) + Sqrt((d*(m + n/2.))/(-m + n/2.))))/
	    (-2*(-m + n/2.)*(m + n/2.) + n) + ((Sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.)) +
	        Sqrt((d*(-m + n/2.))/(m + n/2.)))*((((1 + d - n)*(1 - m - n/2.))/(2.*Power(1 + m - n/2.,2)) -
	           (1 + d - n)/(2.*(1 + m - n/2.)) - (1 - m - n/2.)/(1 + m - n/2.))/
	         (2.*Sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.))) +
	        (d/(2.*(-m + n/2.)) - (d*(m + n/2.))/(2.*Power(-m + n/2.,2)))/(2.*Sqrt((d*(m + n/2.))/(-m + n/2.))))*
	      (-2*d + n))/(-2*(-m + n/2.)*(m + n/2.) + n) +
	   (((-(1 + d - n)/(2.*(1 - m - n/2.)) + ((1 + d - n)*(1 + m - n/2.))/(2.*Power(1 - m - n/2.,2)) -
	           (1 + m - n/2.)/(1 - m - n/2.))/(2.*Sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.))) +
	        (-(d*(-m + n/2.))/(2.*Power(m + n/2.,2)) + d/(2.*(m + n/2.)))/(2.*Sqrt((d*(-m + n/2.))/(m + n/2.))))*
	      (Sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.)) + Sqrt((d*(m + n/2.))/(-m + n/2.)))*(-2*d + n))/
	    (-2*(-m + n/2.)*(m + n/2.) + n);
}

double Func::gamma(double kx, double ky) {
   return 2.*(cos(kx) + cos(ky));
}

double Func::eta(double kx, double ky) {
   return 2.*(cos(kx) - cos(ky));
}

double Func::epsilon(double kx, double ky) {
   return -(gt()*p.t + 3./8.*gs()*p.J*v.chi)*gamma(kx,ky);
}

double Func::BDelta_d() {
   return 3./8.*gs()*p.J*v.BDelta;
}

double Func::BDelta_af() {
   return 2.*gs()*p.J*v.m;
}

double Func::E1(double kx, double ky) {
	double etak = eta(kx, ky);
	double gammak = gamma(kx, ky);
	return -Sqrt(Power(etak,2)*Power((3*BDelta*abgs*J)/8. + 2*lambdaBDelta,2) + Power(lambdam/2. + 2*abgs*J*m,2) +
		     Power(lambdan + mu,2) + Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2) -
		     2*Sqrt(Power(lambdan + mu,2)*(Power(lambdam/2. + 2*abgs*J*m,2) +
		          Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2))));
}

double Func::E3(double kx, double ky) {
	double etak = eta(kx, ky);
	double gammak = gamma(kx, ky);
	return -Sqrt(Power(etak,2)*Power((3*BDelta*abgs*J)/8. + 2*lambdaBDelta,2) + Power(lambdam/2. + 2*abgs*J*m,2) +
		     Power(lambdan + mu,2) + Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2) +
		     2*Sqrt(Power(lambdam/2. + 2*abgs*J*m,2)*Power(lambdan + mu,2) +
		        Power(gammak,2)*Power(lambdan + mu,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2)));
}

double Func::dE1dBDelta(double kx, double ky) {
	double etak = eta(kx, ky);
	double gammak = gamma(kx, ky);
	return  (-3*Power(etak,2)*abgs*J*((3*BDelta*abgs*J)/8. + 2*lambdaBDelta))/
			   (8.*Sqrt(Power(etak,2)*Power((3*BDelta*abgs*J)/8. + 2*lambdaBDelta,2) + Power(lambdam/2. + 2*abgs*J*m,2) +
			       Power(lambdan + mu,2) + Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2) -
			       2*Sqrt(Power(lambdan + mu,2)*(Power(lambdam/2. + 2*abgs*J*m,2) +
			            Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2)))));
}
double Func::dE1dchi(double kx, double ky) {
	double etak = eta(kx, ky);
	double gammak = gamma(kx, ky);
	return -((-3*Power(gammak,2)*abgs*J*((-3*chi*abgs*J)/8. - lambdachi - abgt*t))/4. +
		      (3*Power(gammak,2)*abgs*J*Power(lambdan + mu,2)*((-3*chi*abgs*J)/8. - lambdachi - abgt*t))/
		       (4.*Sqrt(Power(lambdan + mu,2)*(Power(lambdam/2. + 2*abgs*J*m,2) +
		             Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2)))))/
		   (2.*Sqrt(Power(etak,2)*Power((3*BDelta*abgs*J)/8. + 2*lambdaBDelta,2) + Power(lambdam/2. + 2*abgs*J*m,2) +
		       Power(lambdan + mu,2) + Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2) -
		       2*Sqrt(Power(lambdan + mu,2)*(Power(lambdam/2. + 2*abgs*J*m,2) +
		            Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2)))));
}
double Func::dE1dm(double kx, double ky) {
	double etak = eta(kx, ky);
	double gammak = gamma(kx, ky);
	return -(4*abgs*J*(lambdam/2. + 2*abgs*J*m) - (4*abgs*J*(lambdam/2. + 2*abgs*J*m)*Power(lambdan + mu,2))/
		       Sqrt(Power(lambdan + mu,2)*(Power(lambdam/2. + 2*abgs*J*m,2) +
		           Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2))))/
		   (2.*Sqrt(Power(etak,2)*Power((3*BDelta*abgs*J)/8. + 2*lambdaBDelta,2) + Power(lambdam/2. + 2*abgs*J*m,2) +
		       Power(lambdan + mu,2) + Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2) -
		       2*Sqrt(Power(lambdan + mu,2)*(Power(lambdam/2. + 2*abgs*J*m,2) +
		            Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2)))));
}
double Func::dE1dgt(double kx, double ky) {
	double etak = eta(kx, ky);
	double gammak = gamma(kx, ky);
	return -(-2*Power(gammak,2)*t*((-3*chi*abgs*J)/8. - lambdachi - abgt*t) +
		      (2*Power(gammak,2)*Power(lambdan + mu,2)*t*((-3*chi*abgs*J)/8. - lambdachi - abgt*t))/
		       Sqrt(Power(lambdan + mu,2)*(Power(lambdam/2. + 2*abgs*J*m,2) +
		           Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2))))/
		   (2.*Sqrt(Power(etak,2)*Power((3*BDelta*abgs*J)/8. + 2*lambdaBDelta,2) + Power(lambdam/2. + 2*abgs*J*m,2) +
		       Power(lambdan + mu,2) + Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2) -
		       2*Sqrt(Power(lambdan + mu,2)*(Power(lambdam/2. + 2*abgs*J*m,2) +
		            Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2)))));
}
double Func::dE1dgs(double kx, double ky) {
	double etak = eta(kx, ky);
	double gammak = gamma(kx, ky);
	return -((3*BDelta*Power(etak,2)*J*((3*BDelta*abgs*J)/8. + 2*lambdaBDelta))/4. + 4*J*m*(lambdam/2. + 2*abgs*J*m) -
		      (3*chi*Power(gammak,2)*J*((-3*chi*abgs*J)/8. - lambdachi - abgt*t))/4. -
		      (Power(lambdan + mu,2)*(4*J*m*(lambdam/2. + 2*abgs*J*m) -
		           (3*chi*Power(gammak,2)*J*((-3*chi*abgs*J)/8. - lambdachi - abgt*t))/4.))/
		       Sqrt(Power(lambdan + mu,2)*(Power(lambdam/2. + 2*abgs*J*m,2) +
		           Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2))))/
		   (2.*Sqrt(Power(etak,2)*Power((3*BDelta*abgs*J)/8. + 2*lambdaBDelta,2) + Power(lambdam/2. + 2*abgs*J*m,2) +
		       Power(lambdan + mu,2) + Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2) -
		       2*Sqrt(Power(lambdan + mu,2)*(Power(lambdam/2. + 2*abgs*J*m,2) +
		            Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2)))));
}
double Func::dE1dmu(double kx, double ky) {
	double etak = eta(kx, ky);
	double gammak = gamma(kx, ky);
	return -(2*(lambdan + mu) - (2*(lambdan + mu)*(Power(lambdam/2. + 2*abgs*J*m,2) +
	           Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2)))/
	       Sqrt(Power(lambdan + mu,2)*(Power(lambdam/2. + 2*abgs*J*m,2) +
	           Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2))))/
	   (2.*Sqrt(Power(etak,2)*Power((3*BDelta*abgs*J)/8. + 2*lambdaBDelta,2) + Power(lambdam/2. + 2*abgs*J*m,2) +
	       Power(lambdan + mu,2) + Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2) -
	       2*Sqrt(Power(lambdan + mu,2)*(Power(lambdam/2. + 2*abgs*J*m,2) +
	            Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2)))));
}
double Func::dE1dl_chi(double kx, double ky) {
	double etak = eta(kx, ky);
	double gammak = gamma(kx, ky);
	return -(-2*Power(gammak,2)*((-3*chi*abgs*J)/8. - lambdachi - abgt*t) +
		      (2*Power(gammak,2)*Power(lambdan + mu,2)*((-3*chi*abgs*J)/8. - lambdachi - abgt*t))/
		       Sqrt(Power(lambdan + mu,2)*(Power(lambdam/2. + 2*abgs*J*m,2) +
		           Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2))))/
		   (2.*Sqrt(Power(etak,2)*Power((3*BDelta*abgs*J)/8. + 2*lambdaBDelta,2) + Power(lambdam/2. + 2*abgs*J*m,2) +
		       Power(lambdan + mu,2) + Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2) -
		       2*Sqrt(Power(lambdan + mu,2)*(Power(lambdam/2. + 2*abgs*J*m,2) +
		            Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2)))));
}
double Func::dE1dl_n(double kx, double ky) {
	double etak = eta(kx, ky);
	double gammak = gamma(kx, ky);
	return -(2*(lambdan + mu) - (2*(lambdan + mu)*(Power(lambdam/2. + 2*abgs*J*m,2) +
	           Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2)))/
	       Sqrt(Power(lambdan + mu,2)*(Power(lambdam/2. + 2*abgs*J*m,2) +
	           Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2))))/
	   (2.*Sqrt(Power(etak,2)*Power((3*BDelta*abgs*J)/8. + 2*lambdaBDelta,2) + Power(lambdam/2. + 2*abgs*J*m,2) +
	       Power(lambdan + mu,2) + Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2) -
	       2*Sqrt(Power(lambdan + mu,2)*(Power(lambdam/2. + 2*abgs*J*m,2) +
	            Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2)))));
}
double Func::dE1dl_m(double kx, double ky) {
	double etak = eta(kx, ky);
	double gammak = gamma(kx, ky);
	return -(lambdam/2. + 2*abgs*J*m - ((lambdam/2. + 2*abgs*J*m)*Power(lambdan + mu,2))/
		       Sqrt(Power(lambdan + mu,2)*(Power(lambdam/2. + 2*abgs*J*m,2) +
		           Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2))))/
		   (2.*Sqrt(Power(etak,2)*Power((3*BDelta*abgs*J)/8. + 2*lambdaBDelta,2) + Power(lambdam/2. + 2*abgs*J*m,2) +
		       Power(lambdan + mu,2) + Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2) -
		       2*Sqrt(Power(lambdan + mu,2)*(Power(lambdam/2. + 2*abgs*J*m,2) +
		            Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2)))));
}
double Func::dE1dl_BDelta(double kx, double ky) {
	double etak = eta(kx, ky);
	double gammak = gamma(kx, ky);
	return (-2*Power(etak,2)*((3*BDelta*abgs*J)/8. + 2*lambdaBDelta))/
			   Sqrt(Power(etak,2)*Power((3*BDelta*abgs*J)/8. + 2*lambdaBDelta,2) + Power(lambdam/2. + 2*abgs*J*m,2) +
			     Power(lambdan + mu,2) + Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2) -
			     2*Sqrt(Power(lambdan + mu,2)*(Power(lambdam/2. + 2*abgs*J*m,2) +
			          Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2))));
}

double Func::dE3dBDelta(double kx, double ky) {
	double etak = eta(kx, ky);
	double gammak = gamma(kx, ky);
	return (-3*Power(etak,2)*abgs*J*((3*BDelta*abgs*J)/8. + 2*lambdaBDelta))/
			   (8.*Sqrt(Power(etak,2)*Power((3*BDelta*abgs*J)/8. + 2*lambdaBDelta,2) + Power(lambdam/2. + 2*abgs*J*m,2) +
			       Power(lambdan + mu,2) + Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2) +
			       2*Sqrt(Power(lambdam/2. + 2*abgs*J*m,2)*Power(lambdan + mu,2) +
			          Power(gammak,2)*Power(lambdan + mu,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2))));
}
double Func::dE3dchi(double kx, double ky) {
	double etak = eta(kx, ky);
	double gammak = gamma(kx, ky);
	return -((-3*Power(gammak,2)*abgs*J*((-3*chi*abgs*J)/8. - lambdachi - abgt*t))/4. -
		      (3*Power(gammak,2)*abgs*J*Power(lambdan + mu,2)*((-3*chi*abgs*J)/8. - lambdachi - abgt*t))/
		       (4.*Sqrt(Power(lambdam/2. + 2*abgs*J*m,2)*Power(lambdan + mu,2) +
		           Power(gammak,2)*Power(lambdan + mu,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2))))/
		   (2.*Sqrt(Power(etak,2)*Power((3*BDelta*abgs*J)/8. + 2*lambdaBDelta,2) + Power(lambdam/2. + 2*abgs*J*m,2) +
		       Power(lambdan + mu,2) + Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2) +
		       2*Sqrt(Power(lambdam/2. + 2*abgs*J*m,2)*Power(lambdan + mu,2) +
		          Power(gammak,2)*Power(lambdan + mu,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2))));
}
double Func::dE3dm(double kx, double ky) {
	double etak = eta(kx, ky);
	double gammak = gamma(kx, ky);
	return -(4*abgs*J*(lambdam/2. + 2*abgs*J*m) + (4*abgs*J*(lambdam/2. + 2*abgs*J*m)*Power(lambdan + mu,2))/
		       Sqrt(Power(lambdam/2. + 2*abgs*J*m,2)*Power(lambdan + mu,2) +
		         Power(gammak,2)*Power(lambdan + mu,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2)))/
		   (2.*Sqrt(Power(etak,2)*Power((3*BDelta*abgs*J)/8. + 2*lambdaBDelta,2) + Power(lambdam/2. + 2*abgs*J*m,2) +
		       Power(lambdan + mu,2) + Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2) +
		       2*Sqrt(Power(lambdam/2. + 2*abgs*J*m,2)*Power(lambdan + mu,2) +
		          Power(gammak,2)*Power(lambdan + mu,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2))));
}
double Func::dE3dgt(double kx, double ky) {
	double etak = eta(kx, ky);
	double gammak = gamma(kx, ky);
	return -(-2*Power(gammak,2)*t*((-3*chi*abgs*J)/8. - lambdachi - abgt*t) -
		      (2*Power(gammak,2)*Power(lambdan + mu,2)*t*((-3*chi*abgs*J)/8. - lambdachi - abgt*t))/
		       Sqrt(Power(lambdam/2. + 2*abgs*J*m,2)*Power(lambdan + mu,2) +
		         Power(gammak,2)*Power(lambdan + mu,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2)))/
		   (2.*Sqrt(Power(etak,2)*Power((3*BDelta*abgs*J)/8. + 2*lambdaBDelta,2) + Power(lambdam/2. + 2*abgs*J*m,2) +
		       Power(lambdan + mu,2) + Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2) +
		       2*Sqrt(Power(lambdam/2. + 2*abgs*J*m,2)*Power(lambdan + mu,2) +
		          Power(gammak,2)*Power(lambdan + mu,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2))));
}
double Func::dE3dgs(double kx, double ky) {
	double etak = eta(kx, ky);
	double gammak = gamma(kx, ky);
	return -((3*BDelta*Power(etak,2)*J*((3*BDelta*abgs*J)/8. + 2*lambdaBDelta))/4. + 4*J*m*(lambdam/2. + 2*abgs*J*m) -
		      (3*chi*Power(gammak,2)*J*((-3*chi*abgs*J)/8. - lambdachi - abgt*t))/4. +
		      (4*J*m*(lambdam/2. + 2*abgs*J*m)*Power(lambdan + mu,2) -
		         (3*chi*Power(gammak,2)*J*Power(lambdan + mu,2)*((-3*chi*abgs*J)/8. - lambdachi - abgt*t))/4.)/
		       Sqrt(Power(lambdam/2. + 2*abgs*J*m,2)*Power(lambdan + mu,2) +
		         Power(gammak,2)*Power(lambdan + mu,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2)))/
		   (2.*Sqrt(Power(etak,2)*Power((3*BDelta*abgs*J)/8. + 2*lambdaBDelta,2) + Power(lambdam/2. + 2*abgs*J*m,2) +
		       Power(lambdan + mu,2) + Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2) +
		       2*Sqrt(Power(lambdam/2. + 2*abgs*J*m,2)*Power(lambdan + mu,2) +
		          Power(gammak,2)*Power(lambdan + mu,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2))));
}
double Func::dE3dmu(double kx, double ky) {
	double etak = eta(kx, ky);
	double gammak = gamma(kx, ky);
	return -(2*(lambdan + mu) + (2*Power(lambdam/2. + 2*abgs*J*m,2)*(lambdan + mu) +
	         2*Power(gammak,2)*(lambdan + mu)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2))/
	       Sqrt(Power(lambdam/2. + 2*abgs*J*m,2)*Power(lambdan + mu,2) +
	         Power(gammak,2)*Power(lambdan + mu,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2)))/
	   (2.*Sqrt(Power(etak,2)*Power((3*BDelta*abgs*J)/8. + 2*lambdaBDelta,2) + Power(lambdam/2. + 2*abgs*J*m,2) +
	       Power(lambdan + mu,2) + Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2) +
	       2*Sqrt(Power(lambdam/2. + 2*abgs*J*m,2)*Power(lambdan + mu,2) +
	          Power(gammak,2)*Power(lambdan + mu,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2))));
}
double Func::dE3dl_chi(double kx, double ky) {
	double etak = eta(kx, ky);
	double gammak = gamma(kx, ky);
	return -(-2*Power(gammak,2)*((-3*chi*abgs*J)/8. - lambdachi - abgt*t) -
		      (2*Power(gammak,2)*Power(lambdan + mu,2)*((-3*chi*abgs*J)/8. - lambdachi - abgt*t))/
		       Sqrt(Power(lambdam/2. + 2*abgs*J*m,2)*Power(lambdan + mu,2) +
		         Power(gammak,2)*Power(lambdan + mu,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2)))/
		   (2.*Sqrt(Power(etak,2)*Power((3*BDelta*abgs*J)/8. + 2*lambdaBDelta,2) + Power(lambdam/2. + 2*abgs*J*m,2) +
		       Power(lambdan + mu,2) + Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2) +
		       2*Sqrt(Power(lambdam/2. + 2*abgs*J*m,2)*Power(lambdan + mu,2) +
		          Power(gammak,2)*Power(lambdan + mu,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2))));
}
double Func::dE3dl_n(double kx, double ky) {
	double etak = eta(kx, ky);
	double gammak = gamma(kx, ky);
	return -(2*(lambdan + mu) + (2*Power(lambdam/2. + 2*abgs*J*m,2)*(lambdan + mu) +
	         2*Power(gammak,2)*(lambdan + mu)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2))/
	       Sqrt(Power(lambdam/2. + 2*abgs*J*m,2)*Power(lambdan + mu,2) +
	         Power(gammak,2)*Power(lambdan + mu,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2)))/
	   (2.*Sqrt(Power(etak,2)*Power((3*BDelta*abgs*J)/8. + 2*lambdaBDelta,2) + Power(lambdam/2. + 2*abgs*J*m,2) +
	       Power(lambdan + mu,2) + Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2) +
	       2*Sqrt(Power(lambdam/2. + 2*abgs*J*m,2)*Power(lambdan + mu,2) +
	          Power(gammak,2)*Power(lambdan + mu,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2))));
}
double Func::dE3dl_m(double kx, double ky) {
	double etak = eta(kx, ky);
	double gammak = gamma(kx, ky);
	return -(lambdam/2. + 2*abgs*J*m + ((lambdam/2. + 2*abgs*J*m)*Power(lambdan + mu,2))/
		       Sqrt(Power(lambdam/2. + 2*abgs*J*m,2)*Power(lambdan + mu,2) +
		         Power(gammak,2)*Power(lambdan + mu,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2)))/
		   (2.*Sqrt(Power(etak,2)*Power((3*BDelta*abgs*J)/8. + 2*lambdaBDelta,2) + Power(lambdam/2. + 2*abgs*J*m,2) +
		       Power(lambdan + mu,2) + Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2) +
		       2*Sqrt(Power(lambdam/2. + 2*abgs*J*m,2)*Power(lambdan + mu,2) +
		          Power(gammak,2)*Power(lambdan + mu,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2))));
}
double Func::dE3dl_BDelta(double kx, double ky) {
	double etak = eta(kx, ky);
	double gammak = gamma(kx, ky);
	return (-2*Power(etak,2)*((3*BDelta*abgs*J)/8. + 2*lambdaBDelta))/
			   Sqrt(Power(etak,2)*Power((3*BDelta*abgs*J)/8. + 2*lambdaBDelta,2) + Power(lambdam/2. + 2*abgs*J*m,2) +
			     Power(lambdan + mu,2) + Power(gammak,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2) +
			     2*Sqrt(Power(lambdam/2. + 2*abgs*J*m,2)*Power(lambdan + mu,2) +
			        Power(gammak,2)*Power(lambdan + mu,2)*Power((-3*chi*abgs*J)/8. - lambdachi - abgt*t,2)));
}

/*
 *********************************************************
 ******************* CLASS Equations *********************
 *********************************************************
*/

class Equa {
private:
	parameters p; // parameters
	variables v; // variables
	Func *f;

	// result of f (independent on k_x and k_y)
	double gs;
	double dgtdm;
	double dgtdd;
	double dgtdn;
	double dgsdm;
	double dgsdd;
	double dgsdn;
	double BDelta_d;
	double BDelta_af;

	void onePoint(variables &y, int i, int j);

	void refreshRF() {
		gs = f->gs();
		dgtdm = f->dgtdm();
		dgtdd = f->dgtdd();
		dgtdn = f->dgtdn();
		dgsdm = f->dgsdm();
		dgsdd = f->dgsdd();
		dgsdn = f->dgsdn();
		BDelta_d = f->BDelta_d();
		BDelta_af = f->BDelta_af();
	}

public:
	Equa(parameters &p, variables &v) : p(p), v(v) { f = new Func(p, v); refreshRF(); }
	~Equa() { delete f; }

	void setPar(parameters &p) {
		this->p = p;
		f->setPar(p);
		refreshRF();
	}

	void run(variables &w);
};

void Equa::onePoint(variables &w, int i, int j) {
	double x = (double)i*2.*M_PI/p.L;
	double y = (double)j*2.*M_PI/p.L;

	double E[4] = {f->E1(x,y), -f->E1(x,y), f->E3(x,y), -f->E3(x,y) };

	double dEdBDelta[4] = { f->dE1dBDelta(x, y), -f->dE1dBDelta(x, y), f->dE3dBDelta(x, y), -f->dE3dBDelta(x, y) };
	double dEdchi[4] = { f->dE1dchi(x, y), -f->dE1dchi(x, y), f->dE3dchi(x, y), -f->dE3dchi(x, y) };
	double dEdm[4] = { f->dE1dm(x, y), -f->dE1dm(x, y), f->dE3dm(x, y), -f->dE3dm(x, y) };
	double dEdgt[4] = { f->dE1dgt(x, y), -f->dE1dgt(x, y), f->dE3dgt(x, y), -f->dE3dgt(x, y) };
	double dEdgs[4] = { f->dE1dgs(x, y), -f->dE1dgs(x, y), f->dE3dgs(x, y), -f->dE3dgs(x, y) };
	//double dEdmu[4] = { f->dE1dmu(x, y), -f->dE1dmu(x, y), f->dE3dmu(x, y), -f->dE3dmu(x, y) };
	double dEdl_chi[4] = { f->dE1dl_chi(x, y), -f->dE1dl_chi(x, y), f->dE3dl_chi(x, y), -f->dE3dl_chi(x, y)};
	double dEdl_n[4] = { f->dE1dl_n(x, y), -f->dE1dl_n(x, y), f->dE3dl_n(x, y), -f->dE3dl_n(x, y) };
	double dEdl_m[4] = { f->dE1dl_m(x, y), -f->dE1dl_m(x, y), f->dE3dl_m(x, y), -f->dE3dl_m(x, y) };
	double dEdl_BDelta[4] = { f->dE1dl_BDelta(x, y), -f->dE1dl_BDelta(x, y), f->dE3dl_BDelta(x, y), -f->dE3dl_BDelta(x, y) };

	double Etan[4] = { (1-tanh(p.beta*E[0]/2.))/2., (1-tanh(p.beta*E[1]/2.))/2., (1-tanh(p.beta*E[2]/2.))/2., (1-tanh(p.beta*E[3]/2.))/2. };

	for(int i=0; i<4; i++) {
		w.BDelta += dEdBDelta[i]*Etan[i];
		w.chi += dEdchi[i]*Etan[i];
		w.m += (dEdm[i] + dEdgt[i]*dgtdm + dEdgs[i]*dgsdm)*Etan[i];
		w.d += (dEdgt[i]*dgtdd + dEdgs[i]*dgsdd)*Etan[i];
		w.mu += (-dEdgt[i]*dgtdn - dEdgs[i]*dgsdn)*Etan[i];
		w.l_BDelta += dEdl_BDelta[i]*Etan[i];
		w.l_chi += dEdl_chi[i]*Etan[i];
		w.l_m += dEdl_m[i]*Etan[i];
		w.l_n += dEdl_n[i]*Etan[i];
	}
}

void Equa::run(variables &w)
{
	clear(w);
	variables temp; clear(temp); // temporary vector using during summation
	unsigned int max = p.L/2; // attention: L should be even!

//	int nthreads=1;
//	#pragma omp parallel
//	{
//		nthreads = omp_get_num_threads();
//	}

	// SERIAL CASE
//	if(nthreads != 8) {
        // SMART WAY OF SUMMATION - I use symmetry k_x,y <-> -k_x,y
        // * edges (2 times)
        for(unsigned int i=1; i<max; i++) {
    		onePoint(temp, i, 0);
    		onePoint(temp, 0, i);
    		onePoint(temp, i, max-i);
        }
        addTo(temp, w, 2);
        clear(temp);
        // * interior area (4 times) - SERIAL CASE
        for(unsigned int i=1; i<max; i++) {
        	for(unsigned int j=1; j<i; j++) {
        		onePoint(temp, max-i, j);
        	}
        }
        addTo(temp, w, 4);
        clear(temp);
        // * missing corners
        onePoint(w, max, 0);
        onePoint(w, 0, max);
        onePoint(w, 0, 0);
//	}
//	else { // PARALEL CASE ONLY IF 8 THREADS: export OMP_NUM_THREADS=8
//		if(max%8 != 0) { fprintf(stderr, "max mod 8 != 0\n"); exit(EXIT_FAILURE); }
//		// SMART WAY OF SUMMATION - I use symmetry k_x,y <-> -k_x,y
//		// * edges (2 times)
//		for(unsigned int i=1; i<max; i++) {
//			onePoint(temp, i, 0);
//			onePoint(temp, 0, i);
//			onePoint(temp, i, max-i);
//		}
//		addTo(temp, w, 2.);
//		clear(temp);
//
//		// * interior area (4 times)
//		unsigned int th_id, th_start, th_end;
//		variables th_temp[8];
//		for(int i=0; i<8; i++) { clear(th_temp[i]); }
//		#pragma omp parallel private(th_id, th_start, th_end)
//		{
//			th_id = omp_get_thread_num();
//			th_start = (unsigned)(((double)max)*(4.-sqrt(16.-2.*th_id))/4.);
//			th_end = (unsigned)(((double)max)*(4.-sqrt(16.-2.*(th_id+1.)))/4.);
//			for(unsigned int i=th_start; i<th_end; i++) {
//				for(unsigned int j=1; j<i; j++) {
//			           onePoint(th_temp[th_id], max-i, j);
//				}
//			}
//		}
//		for(int i=0; i<8; i++) {
//			temp.BDelta += th_temp[i].BDelta;
//			temp.chi  += th_temp[i].chi;
//			temp.m  += th_temp[i].m;
//			temp.d  += th_temp[i].d;
//			temp.mu  += th_temp[i].mu;
//			temp.l_BDelta  += th_temp[i].l_BDelta;
//			temp.l_chi  += th_temp[i].l_chi;
//			temp.l_m  += th_temp[i].l_m;
//			temp.l_n  += th_temp[i].l_n;
//		}
//		addTo(temp, w, 4);
//		clear(temp);
//		// * missing corners
//		onePoint(w, max, 0);
//		onePoint(w, 0, max);
//		onePoint(w, 0, 0);
//	}

	double L2 = (double)(p.L*p.L);
	w.BDelta /= L2;
	w.chi /= L2;
	w.m /= L2;
	w.d /= L2;
	w.mu /= L2;
	w.l_BDelta /= L2;
	w.l_chi /= L2;
	w.l_m /= L2;
	w.l_n /= L2;

	w.BDelta += 4*v.l_BDelta + 3./2.*gs*p.J*v.BDelta;
	w.chi += 4*v.l_chi + 3./2.*gs*p.J*v.chi;
	w.m += v.l_m + 4*gs*p.J*v.m;
	w.d += p.U;
	w.mu += -v.l_n; //-v.mu - v.l_n;
	w.l_BDelta += 4*v.BDelta;
	w.l_chi += 4*v.chi;
	w.l_m += v.m;
	w.l_n += -p.sdelta;

	//if(v.m < 0) w.m = v.m-0.4; // temporary for debuging
	if(w.BDelta != w.BDelta) { w.BDelta = 1000; }
	if(w.chi != w.chi) { w.chi = 1000; }
	if(w.m != w.m) { w.m = 1000; }
	if(w.d != w.d) { w.d = -10*v.d; }
	if(w.mu != w.mu) { w.mu = 1000; }
	if(w.l_BDelta != w.l_BDelta) { w.l_BDelta = 1000; }
	if(w.l_chi != w.l_chi) { w.l_chi = 1000; }
	if(w.l_m != w.l_m) { w.l_m = 1000; }
	if(w.l_n != w.l_n) { w.l_n = 1000; }
}

